Spatial regression analysis for the total number of cases

load("../../data/unionedListw_d_berlinNotSubdivided.Rdata")
load("../../data/kreiseWithCovidMeldedatumWeeklyPredictors.Rdata")

Explorative analysis

Scatterplot matrix

st_drop_geometry(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded) %>%
    dplyr::select(incidenceTotal, Langzeitarbeitslosenquote, Hochqualifizierte, Breitbandversorgung,
        Stimmenanteile.AfD, Studierende, Einpendler, Auspendler, Einpendler, Auslaenderanteil,
        Laendlichkeit) %>%
    ggpairs()

We see some predictors correlating with the incidence rate as well as some moderate collinearity between some of our predictors.

Maps of response and predictors

For reference I create a series of maps for the response and all predictors.

kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
    tm_shape() + tm_polygons(col = c("incidenceTotal"), legend.hist = TRUE, palette = "-plasma",
    legend.reverse = TRUE, title = "Incidence rate for total period") + tm_layout(legend.outside = TRUE,
    bg.color = "darkgrey", outer.bg.color = "lightgrey", legend.outside.size = 0.5,
    attr.outside = TRUE, main.title = "Incidence rate total period") + tm_scale_bar()

mForeigner <- kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
    tm_shape() + tm_polygons(col = c("Auslaenderanteil"), legend.hist = TRUE, legend.reverse = TRUE,
    title = "Share of foreigners", style = "pretty") + tm_layout(legend.outside = TRUE,
    bg.color = "darkgrey", outer.bg.color = "lightgrey", legend.outside.size = 0.5,
    attr.outside = TRUE, main.title = "Foreign population") + tm_scale_bar()

print(mForeigner)

kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
    tm_shape() + tm_polygons(col = c("Auspendler"), legend.hist = TRUE, legend.reverse = TRUE,
    title = "Outgoing commuters as share\nof total employees", style = "kmeans") +
    tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
        legend.outside.size = 0.5, attr.outside = TRUE, main.title = "Outgoing commuters") +
    tm_scale_bar()

kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
    tm_shape() + tm_polygons(col = c("Einpendler"), legend.hist = TRUE, legend.reverse = TRUE,
    title = "Incomming commuters as share\nof total employees", style = "kmeans") +
    tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
        legend.outside.size = 0.5, attr.outside = TRUE, main.title = "Incomming commuters") +
    tm_scale_bar()

mRural <- kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
    tm_shape() + tm_polygons(col = c("Laendlichkeit"), legend.hist = TRUE, palette = "Purples",
    legend.reverse = TRUE, title = "Share of inhabitants in places\nwith less then 150 Inh/sqkm",
    style = "kmeans") + tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
    legend.outside.size = 0.5, attr.outside = TRUE, main.title = "Ruralness") + tm_scale_bar()
print(mRural)

kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
    tm_shape() + tm_polygons(col = c("Studierende"), legend.hist = TRUE, palette = "Purples",
    legend.reverse = TRUE, title = "Students per \n1000 Inhabitants", style = "kmeans") +
    tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
        legend.outside.size = 0.5, attr.outside = TRUE, main.title = "Students") +
    tm_scale_bar()

kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
    tm_shape() + tm_polygons(col = c("Hochqualifizierte"), legend.hist = TRUE, legend.reverse = TRUE,
    title = "Highly qualified employees/ total employees", style = "kmeans") + tm_layout(legend.outside = TRUE,
    bg.color = "darkgrey", outer.bg.color = "lightgrey", legend.outside.size = 0.5,
    attr.outside = TRUE, main.title = "Highly qualified employees") + tm_scale_bar()

kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
    tm_shape() + tm_polygons(col = c("Breitbandversorgung"), legend.hist = TRUE,
    palette = "Purples", legend.reverse = TRUE, title = "Share of households with internet \nconnectivity > 5 mBits/s",
    style = "kmeans") + tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
    legend.outside.size = 0.5, attr.outside = TRUE, main.title = "Highspeed internet") +
    tm_scale_bar()

kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
    tm_shape() + tm_polygons(col = c("Langzeitarbeitslosenquote"), legend.hist = TRUE,
    palette = "Oranges", legend.reverse = TRUE, title = "Unemployment rate [%]") +
    tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
        legend.outside.size = 0.5, attr.outside = TRUE, main.title = "Long term unemloyment rate") +
    tm_scale_bar()

kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
    tm_shape() + tm_polygons(col = c("Stimmenanteile.AfD"), legend.hist = TRUE, palette = "Blues",
    legend.reverse = TRUE, title = "Share of votes") + tm_layout(legend.outside = TRUE,
    bg.color = "darkgrey", outer.bg.color = "lightgrey", legend.outside.size = 0.5,
    attr.outside = TRUE, main.title = "Right winged party (AfD)") + tm_scale_bar()

Normal GLM

We start with a normal GLM before checking for spatial autocorrelation in the residuals. Since we have count data a Poisson GLM with an offset for the population at risk seems a natural choice.

Poisson GLM

modelGlm <- glm(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
    Langzeitarbeitslosenquote + Einpendler + Studierende + Laendlichkeit + offset(log(EWZ)),
    data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded, family = poisson)
summary(modelGlm)

Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + 
    Hochqualifizierte + Langzeitarbeitslosenquote + Einpendler + 
    Studierende + Laendlichkeit + offset(log(EWZ)), family = poisson, 
    data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-75.552  -13.633   -2.309   10.382   70.683  

Coefficients:
                            Estimate Std. Error z value Pr(>|z|)    
(Intercept)               -3.449e+00  4.460e-03 -773.31   <2e-16 ***
Stimmenanteile.AfD         4.122e-02  9.852e-05  418.38   <2e-16 ***
Auslaenderanteil           3.289e-02  1.226e-04  268.35   <2e-16 ***
Hochqualifizierte         -1.276e-02  1.185e-04 -107.66   <2e-16 ***
Langzeitarbeitslosenquote -2.364e-02  3.704e-04  -63.81   <2e-16 ***
Einpendler                -1.599e-03  4.323e-05  -36.99   <2e-16 ***
Studierende                1.509e-04  1.211e-05   12.46   <2e-16 ***
Laendlichkeit             -1.327e-03  2.565e-05  -51.72   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 344479  on 399  degrees of freedom
Residual deviance: 156226  on 392  degrees of freedom
AIC: 160590

Number of Fisher Scoring iterations: 4

Negative binomial GLM to account for overdispersion

Since we should be suspicious with respect to overdispersion we will run a negative binomial and afterwards a quasi-poisson GLM to account for that. Since the negative binomial GLM triggers some complications when using it with the spatial eigenvector mapping we will stay with the quasi-poisson model afterwards.

modelGlmNb <- glm.nb(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
    Langzeitarbeitslosenquote + Einpendler + Studierende + Laendlichkeit + offset(log(EWZ)),
    data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded)
summary(modelGlmNb)

Call:
glm.nb(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + 
    Hochqualifizierte + Langzeitarbeitslosenquote + Einpendler + 
    Studierende + Laendlichkeit + offset(log(EWZ)), data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded, 
    init.theta = 21.58561699, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.4933  -0.7228  -0.0782   0.5668   3.1286  

Coefficients:
                            Estimate Std. Error z value Pr(>|z|)    
(Intercept)               -3.5407641  0.1201151 -29.478  < 2e-16 ***
Stimmenanteile.AfD         0.0458598  0.0022967  19.968  < 2e-16 ***
Auslaenderanteil           0.0376844  0.0030124  12.510  < 2e-16 ***
Hochqualifizierte         -0.0144682  0.0029854  -4.846 1.26e-06 ***
Langzeitarbeitslosenquote -0.0480101  0.0090316  -5.316 1.06e-07 ***
Einpendler                -0.0012145  0.0013595  -0.893    0.372    
Studierende                0.0002775  0.0002649   1.048    0.295    
Laendlichkeit             -0.0008254  0.0005037  -1.639    0.101    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(21.5856) family taken to be 1)

    Null deviance: 851.01  on 399  degrees of freedom
Residual deviance: 403.29  on 392  degrees of freedom
AIC: 7156.5

Number of Fisher Scoring iterations: 1

              Theta:  21.59 
          Std. Err.:  1.52 

 2 x log-likelihood:  -7138.506 
drop1(modelGlmNb)
Single term deletions

Model:
sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte + 
    Langzeitarbeitslosenquote + Einpendler + Studierende + Laendlichkeit + 
    offset(log(EWZ))
                          Df Deviance    AIC
<none>                         403.29 7154.5
Stimmenanteile.AfD         1   812.34 7561.6
Auslaenderanteil           1   556.93 7306.1
Hochqualifizierte          1   427.05 7176.3
Langzeitarbeitslosenquote  1   431.45 7180.7
Einpendler                 1   404.07 7153.3
Studierende                1   404.45 7153.7
Laendlichkeit              1   405.90 7155.1
modelGlmNb <- update(modelGlmNb, ~. - Einpendler)
summary(modelGlmNb)

Call:
glm.nb(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + 
    Hochqualifizierte + Langzeitarbeitslosenquote + Studierende + 
    Laendlichkeit + offset(log(EWZ)), data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded, 
    init.theta = 21.54407359, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.4684  -0.7223  -0.1059   0.5948   3.1251  

Coefficients:
                            Estimate Std. Error z value Pr(>|z|)    
(Intercept)               -3.6325306  0.0626247 -58.005  < 2e-16 ***
Stimmenanteile.AfD         0.0458152  0.0022989  19.929  < 2e-16 ***
Auslaenderanteil           0.0377287  0.0030153  12.513  < 2e-16 ***
Hochqualifizierte         -0.0142960  0.0029816  -4.795 1.63e-06 ***
Langzeitarbeitslosenquote -0.0437307  0.0076450  -5.720 1.06e-08 ***
Studierende                0.0003129  0.0002617   1.196    0.232    
Laendlichkeit             -0.0008165  0.0005042  -1.619    0.105    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(21.5441) family taken to be 1)

    Null deviance: 849.38  on 399  degrees of freedom
Residual deviance: 403.30  on 393  degrees of freedom
AIC: 7155.3

Number of Fisher Scoring iterations: 1

              Theta:  21.54 
          Std. Err.:  1.52 

 2 x log-likelihood:  -7139.284 
modelGlmNb <- update(modelGlmNb, ~. - Studierende)
summary(modelGlmNb)

Call:
glm.nb(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + 
    Hochqualifizierte + Langzeitarbeitslosenquote + Laendlichkeit + 
    offset(log(EWZ)), data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded, 
    init.theta = 21.46423243, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.4607  -0.7231  -0.0854   0.5758   3.1371  

Coefficients:
                            Estimate Std. Error z value Pr(>|z|)    
(Intercept)               -3.6403759  0.0624123 -58.328  < 2e-16 ***
Stimmenanteile.AfD         0.0454243  0.0022782  19.939  < 2e-16 ***
Auslaenderanteil           0.0377774  0.0030208  12.506  < 2e-16 ***
Hochqualifizierte         -0.0126774  0.0026518  -4.781 1.75e-06 ***
Langzeitarbeitslosenquote -0.0424379  0.0075784  -5.600 2.15e-08 ***
Laendlichkeit             -0.0008477  0.0005042  -1.681   0.0927 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(21.4642) family taken to be 1)

    Null deviance: 846.24  on 399  degrees of freedom
Residual deviance: 403.31  on 394  degrees of freedom
AIC: 7154.8

Number of Fisher Scoring iterations: 1

              Theta:  21.46 
          Std. Err.:  1.51 

 2 x log-likelihood:  -7140.788 

Quasi-Poisson GLM to account for overdispersion

modelGlmQp <- glm(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
    Langzeitarbeitslosenquote + Einpendler + Studierende + Laendlichkeit + offset(log(EWZ)),
    data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded, family = quasipoisson)
summary(modelGlmQp)

Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + 
    Hochqualifizierte + Langzeitarbeitslosenquote + Einpendler + 
    Studierende + Laendlichkeit + offset(log(EWZ)), family = quasipoisson, 
    data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-75.552  -13.633   -2.309   10.382   70.683  

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)               -3.4489644  0.0884963 -38.973  < 2e-16 ***
Stimmenanteile.AfD         0.0412170  0.0019548  21.085  < 2e-16 ***
Auslaenderanteil           0.0328903  0.0024319  13.524  < 2e-16 ***
Hochqualifizierte         -0.0127585  0.0023516  -5.426 1.01e-07 ***
Langzeitarbeitslosenquote -0.0236376  0.0073502  -3.216  0.00141 ** 
Einpendler                -0.0015993  0.0008578  -1.864  0.06300 .  
Studierende                0.0001509  0.0002402   0.628  0.53040    
Laendlichkeit             -0.0013267  0.0005090  -2.607  0.00949 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 393.7124)

    Null deviance: 344479  on 399  degrees of freedom
Residual deviance: 156226  on 392  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4
drop1(modelGlmQp, test = "F")
Single term deletions

Model:
sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte + 
    Langzeitarbeitslosenquote + Einpendler + Studierende + Laendlichkeit + 
    offset(log(EWZ))
                          Df Deviance  F value    Pr(>F)    
<none>                         156226                       
Stimmenanteile.AfD         1   322415 416.9966 < 2.2e-16 ***
Auslaenderanteil           1   225877 174.7648 < 2.2e-16 ***
Hochqualifizierte          1   167833  29.1225 1.179e-07 ***
Langzeitarbeitslosenquote  1   160333  10.3051  0.001436 ** 
Einpendler                 1   157593   3.4279  0.064854 .  
Studierende                1   156381   0.3869  0.534276    
Laendlichkeit              1   158920   6.7577  0.009687 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelGlmQp <- update(modelGlmQp, ~. - Studierende)
drop1(modelGlmQp, test = "F")
Single term deletions

Model:
sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte + 
    Langzeitarbeitslosenquote + Einpendler + Laendlichkeit + 
    offset(log(EWZ))
                          Df Deviance  F value    Pr(>F)    
<none>                         156381                       
Stimmenanteile.AfD         1   325141 424.1105 < 2.2e-16 ***
Auslaenderanteil           1   225948 174.8296 < 2.2e-16 ***
Hochqualifizierte          1   168755  31.0982  4.57e-08 ***
Langzeitarbeitslosenquote  1   160334   9.9362  0.001745 ** 
Einpendler                 1   157783   3.5239  0.061231 .  
Laendlichkeit              1   159117   6.8761  0.009075 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelGlmQp <- update(modelGlmQp, ~. - Langzeitarbeitslosenquote)
drop1(modelGlmQp, test = "F")
Single term deletions

Model:
sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte + 
    Einpendler + Laendlichkeit + offset(log(EWZ))
                   Df Deviance  F value    Pr(>F)    
<none>                  160334                       
Stimmenanteile.AfD  1   325149 405.0093 < 2.2e-16 ***
Auslaenderanteil    1   229927 171.0135 < 2.2e-16 ***
Hochqualifizierte   1   170097  23.9891 1.417e-06 ***
Einpendler          1   160405   0.1741   0.67669    
Laendlichkeit       1   161683   3.3146   0.06943 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelGlmQp <- update(modelGlmQp, ~. - Einpendler)
drop1(modelGlmQp, test = "F")
Single term deletions

Model:
sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte + 
    Laendlichkeit + offset(log(EWZ))
                   Df Deviance  F value    Pr(>F)    
<none>                  160405                       
Stimmenanteile.AfD  1   325303 406.0623 < 2.2e-16 ***
Auslaenderanteil    1   229951 171.2579 < 2.2e-16 ***
Hochqualifizierte   1   170969  26.0133 5.274e-07 ***
Laendlichkeit       1   161811   3.4608   0.06358 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Explained deviance:

1 - modelGlmNb$deviance/modelGlmNb$null.deviance
[1] 0.5234146
summary(modelGlmQp)

Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + 
    Hochqualifizierte + Laendlichkeit + offset(log(EWZ)), family = quasipoisson, 
    data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-76.612  -13.848   -2.166   11.822   72.191  

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -3.6310778  0.0486591 -74.623  < 2e-16 ***
Stimmenanteile.AfD  0.0399200  0.0019120  20.879  < 2e-16 ***
Auslaenderanteil    0.0327747  0.0024445  13.408  < 2e-16 ***
Hochqualifizierte  -0.0101702  0.0019856  -5.122 4.74e-07 ***
Laendlichkeit      -0.0009186  0.0004915  -1.869   0.0624 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 400.4526)

    Null deviance: 344479  on 399  degrees of freedom
Residual deviance: 160405  on 395  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

Explained deviance:

1 - modelGlmQp$deviance/modelGlmQp$null.deviance
[1] 0.5343536

We end up with a model of decent quality. Directions of effect seem to be aligned with expectations:

  • higher share of votes for right winged party is associated with higher incidence. Presumably a proxy for the share of population opposing mask wearing and social distancing and vaccination
  • higher share of foreigners is associated with higher incidence. Foreigners might not be reach by information campaigns with respect to social distancing due to language problems
  • higher share of highly qualified work force is associated with lower incidence rates. For those employees it might be easier to work from home office and to avoid close contacts during work hours at office
  • higher ruralness (higher share of population living in rural areas) is associated with lower incidence rates. This might be due to lower contact rates e.g. by lower share of public transport.

Checking for spatial autocorrelation in the residuals

Global Moran’s I

For regression residuals we need to use a different test as residuals are centered around zero and will sum up to zero.

lm.morantest(modelGlmNb, listw = unionedListw_d)

    Global Moran I for regression residuals

data:  
model: glm.nb(formula = sumCasesTotal ~ Stimmenanteile.AfD +
Auslaenderanteil + Hochqualifizierte + Langzeitarbeitslosenquote +
Laendlichkeit + offset(log(EWZ)), data =
kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded, init.theta =
21.46423243, link = log)
weights: unionedListw_d

Moran I statistic standard deviate = 22.666, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
    0.5078558227    -0.0003562342     0.0005027223 
(moranGlmQp <- lm.morantest(modelGlmQp, listw = unionedListw_d))

    Global Moran I for regression residuals

data:  
model: glm(formula = sumCasesTotal ~ Stimmenanteile.AfD +
Auslaenderanteil + Hochqualifizierte + Laendlichkeit +
offset(log(EWZ)), family = quasipoisson, data =
kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded)
weights: unionedListw_d

Moran I statistic standard deviate = 22.628, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
    5.079319e-01    -5.494062e-07     5.038494e-04 

Both model suffer from significant global spatial autocorrelation.

This implies that the usual assumption about independence of errors is violated. In turn, our standard errors might be too low, p-values too small, size (and potentially even sign) of the regression coefficients might be wrong. So we need to incorporate spatial autocorrelation in our analysis.

Plot residuals

kreiseCentroids <- st_centroid(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
    of_largest_polygon = TRUE)

Add residuals to centroids for plotting

kreiseCentroids$residGlmNb <- residuals(modelGlmNb)
kreiseCentroids$residGlmNbAbs <- abs(kreiseCentroids$residGlmNb)

kreiseCentroids$residGlmQp <- residuals(modelGlmQp)
kreiseCentroids$residGlmQpAbs <- abs(kreiseCentroids$residGlmQp)
m1 <- tm_shape(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded) + tm_polygons(col = "grey") +
    tm_shape(kreiseCentroids) + tm_bubbles(size = "residGlmNbAbs", col = "residGlmNb",
    palette = "-RdBu", alpha = 0.9, perceptual = TRUE, scale = 0.8, border.alpha = 0.3,
    title.size = "Abs residual", title.col = "Residuals", n = 3) + tm_layout(main.title = "Pearson residuals, GLM NB",
    bg = "darkgrey", legend.outside = TRUE, attr.outside = TRUE) + tm_scale_bar()

m2 <- tm_shape(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded) + tm_polygons(col = "grey") +
    tm_shape(kreiseCentroids) + tm_bubbles(size = "residGlmNbAbs", col = "residGlmNb",
    palette = "-RdBu", alpha = 0.9, perceptual = TRUE, scale = 0.8, border.alpha = 0.3,
    title.size = "Abs residual", title.col = "Residuals", n = 3) + tm_layout(main.title = "Pearson residuals, GLM QP",
    bg = "darkgrey", legend.outside = TRUE, attr.outside = TRUE) + tm_scale_bar()

tmap_arrange(m1, m2)

As indicated by global Moran’s I we see that large positive and large negative residuals form some cluster. The differences between both GLM models seem neglectable.

Spatial Eigenvector Mapping

The idea behind the spatial eigenvector mapping approach is to use additional covariates that aborb the spatial autocorrelation, leading to unbiased estimators for the other predictors. The additional covariates are based on the eigenfunction decomposition of the spatial weight matrix \(W\). Eigenvectors of \(W\) represent the decompositions the spatial weight Matrix into all mutually orthogonal eigenvectors. Those with positive eigenvalues represent positive autocorrelation, whereas eigenvectors with negative eigenvalues represent negative autocorrelation. Only eigenvectors with positive eigenvalues are used for the selection.

Selection of eigenvectors

The function ME uses brute force eigenvector selection to reach a subset of such vectors to be added to the RHS of the GLM model to reduce residual autocorrelation to below the specified alpha value (defaults to 0.05). Since eigenvector selection only works on symmetric weights, the weights are made symmetric beforehand.

meQp <- spatialreg::ME(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
    Laendlichkeit, family = quasipoisson, data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
    offset = log(EWZ), listw = unionedListw_d)

Refitting GLM under incorporation of the selected spatial eigenvectors:

modelSevmQp <- glm(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
    Laendlichkeit + fitted(meQp), family = quasipoisson, data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
    offset = log(EWZ))
summary(modelSevmQp)

Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + 
    Hochqualifizierte + Laendlichkeit + fitted(meQp), family = quasipoisson, 
    data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded, 
    offset = log(EWZ))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-52.528  -10.397   -1.632    6.793   56.002  

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -3.6821280  0.0447482 -82.285  < 2e-16 ***
Stimmenanteile.AfD  0.0405874  0.0018294  22.186  < 2e-16 ***
Auslaenderanteil    0.0286188  0.0022504  12.717  < 2e-16 ***
Hochqualifizierte  -0.0051287  0.0018140  -2.827 0.004941 ** 
Laendlichkeit      -0.0004370  0.0004085  -1.070 0.285383    
fitted(meQp)vec2    1.2129876  0.1909764   6.352 6.03e-10 ***
fitted(meQp)vec6   -0.8132039  0.1744241  -4.662 4.32e-06 ***
fitted(meQp)vec4   -1.7032836  0.1475831 -11.541  < 2e-16 ***
fitted(meQp)vec10   0.1101595  0.1610011   0.684 0.494251    
fitted(meQp)vec1   -0.7785702  0.1576694  -4.938 1.18e-06 ***
fitted(meQp)vec18  -0.7198156  0.1481683  -4.858 1.73e-06 ***
fitted(meQp)vec21  -1.4354292  0.1692912  -8.479 4.95e-16 ***
fitted(meQp)vec11  -0.3051171  0.1477679  -2.065 0.039609 *  
fitted(meQp)vec19  -1.0557902  0.1680532  -6.282 9.04e-10 ***
fitted(meQp)vec5    0.5968183  0.1644262   3.630 0.000322 ***
fitted(meQp)vec8    0.5483723  0.1670307   3.283 0.001121 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 224.0348)

    Null deviance: 344479  on 399  degrees of freedom
Residual deviance:  85686  on 384  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

The procedure added 11 spatial eigenvectors to the model. Together this leads to a more than satisfying amount of explained deviance. However, we need to keep in mind that a good share of that come from the spatial eigenvectors.

1 - modelSevmQp$deviance/modelSevmQp$null.deviance
[1] 0.7512578

Ruralness of the district now became insignificant so we might want to drop it from the model.

meQp <- spatialreg::ME(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte,
    family = quasipoisson, data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
    offset = log(EWZ), listw = unionedListw_d)

Refitting GLM under incorporation of the selected spatial eigenvectors:

modelSevmQp <- glm(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
    fitted(meQp), family = quasipoisson, offset = log(EWZ), data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded)
summary(modelSevmQp)

Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + 
    Hochqualifizierte + fitted(meQp), family = quasipoisson, 
    data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded, 
    offset = log(EWZ))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-53.356  -10.321   -1.947    6.936   56.477  

Coefficients:
                    Estimate Std. Error  t value Pr(>|t|)    
(Intercept)        -3.721873   0.036744 -101.293  < 2e-16 ***
Stimmenanteile.AfD  0.041288   0.001849   22.333  < 2e-16 ***
Auslaenderanteil    0.030151   0.002176   13.853  < 2e-16 ***
Hochqualifizierte  -0.004930   0.001718   -2.870  0.00433 ** 
fitted(meQp)vec2    1.229067   0.193236    6.360 5.69e-10 ***
fitted(meQp)vec6   -0.791620   0.177978   -4.448 1.14e-05 ***
fitted(meQp)vec4   -1.718902   0.149249  -11.517  < 2e-16 ***
fitted(meQp)vec10   0.132188   0.163657    0.808  0.41975    
fitted(meQp)vec1   -0.661805   0.146610   -4.514 8.46e-06 ***
fitted(meQp)vec18  -0.713091   0.149129   -4.782 2.48e-06 ***
fitted(meQp)vec21  -1.467720   0.171605   -8.553 2.85e-16 ***
fitted(meQp)vec11  -0.332178   0.148986   -2.230  0.02635 *  
fitted(meQp)vec19  -1.001783   0.169929   -5.895 8.17e-09 ***
fitted(meQp)vec5    0.549126   0.167337    3.282  0.00113 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 229.6598)

    Null deviance: 344479  on 399  degrees of freedom
Residual deviance:  88520  on 386  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

The eigenvectors selected did not change by excluding ruralness from the model.

Plotting selected spatial eigenvectors

Presumably we have missed a lot of other predictors that are now partially absorbed into different eigenvectors. It might be worth to plot and investigate the eigenvectors that made it into the model. Therefore, we attach the selected eigenvectors to the sf object and plot them.

summary(fitted(meQp))
      vec2               vec6                vec4               vec10          
 Min.   :-0.17702   Min.   :-0.190177   Min.   :-0.107585   Min.   :-0.148895  
 1st Qu.:-0.04621   1st Qu.:-0.014842   1st Qu.:-0.030753   1st Qu.:-0.018094  
 Median : 0.02081   Median : 0.003256   Median :-0.008749   Median : 0.003103  
 Mean   : 0.00000   Mean   : 0.000000   Mean   : 0.000000   Mean   : 0.000000  
 3rd Qu.: 0.03945   3rd Qu.: 0.013008   3rd Qu.: 0.040463   3rd Qu.: 0.017146  
 Max.   : 0.07715   Max.   : 0.185980   Max.   : 0.112574   Max.   : 0.141635  
      vec1               vec18                vec21           
 Min.   :-0.157693   Min.   :-0.2348596   Min.   :-0.1337760  
 1st Qu.:-0.028969   1st Qu.:-0.0206895   1st Qu.:-0.0257369  
 Median : 0.001129   Median : 0.0003868   Median :-0.0003304  
 Mean   : 0.000000   Mean   : 0.0000000   Mean   : 0.0000000  
 3rd Qu.: 0.031472   3rd Qu.: 0.0263951   3rd Qu.: 0.0304928  
 Max.   : 0.115554   Max.   : 0.2985828   Max.   : 0.2642922  
     vec11                vec19                 vec5          
 Min.   :-0.2190326   Min.   :-0.1735843   Min.   :-0.224152  
 1st Qu.:-0.0184901   1st Qu.:-0.0303940   1st Qu.:-0.026723  
 Median : 0.0007727   Median :-0.0002972   Median :-0.007429  
 Mean   : 0.0000000   Mean   : 0.0000000   Mean   : 0.000000  
 3rd Qu.: 0.0151115   3rd Qu.: 0.0299778   3rd Qu.: 0.021318  
 Max.   : 0.1915285   Max.   : 0.1358331   Max.   : 0.162426  
sevQp <- fitted(meQp)
kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem <- st_sf(data.frame(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
    sevQp))
tm_shape(kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem) + tm_polygons(col = colnames(sevQp),
    palette = "-RdBu", lwd = 0.5, n = 6, midpoint = 0, legend.show = FALSE) + tm_layout(main.title = "Selected spatial eigenvectors",
    legend.outside = TRUE, attr.outside = TRUE, panel.show = TRUE, panel.labels = colnames(sevQp)) +
    tm_scale_bar()

Rechecking spatial autocorrelation

(moranSevmQp <- lm.morantest(modelSevmQp, listw = unionedListw_d))

    Global Moran I for regression residuals

data:  
model: glm(formula = sumCasesTotal ~ Stimmenanteile.AfD +
Auslaenderanteil + Hochqualifizierte + fitted(meQp), family =
quasipoisson, data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
offset = log(EWZ))
weights: unionedListw_d

Moran I statistic standard deviate = 7.3156, p-value = 1.281e-13
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
    1.680228e-01    -2.875261e-06     5.275358e-04 

We see that were is still some left over spatial autocorrelation not absorbed by the spatial eigenvectors. However, the amount of spatial autocorrelation has been reduced by a strong degree, from 0.51 to 0.17.

kreiseCentroids$residSevmQp <- residuals(modelSevmQp)
kreiseCentroids$residSevmQpAbs <- abs(kreiseCentroids$residSevmQp)

tm_shape(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded) + tm_polygons(col = "grey") +
    tm_shape(kreiseCentroids) + tm_bubbles(size = "residSevmQpAbs", col = "residSevmQp",
    palette = "-RdBu", alpha = 0.9, perceptual = TRUE, scale = 0.8, border.alpha = 0.3,
    title.size = "Abs residual", title.col = "Residuals", n = 5) + tm_layout(main.title = "Pearson residuals, SEVM NB",
    bg = "darkgrey", legend.outside = TRUE, attr.outside = TRUE) + tm_scale_bar()

## Spatial varying coefficients?

Which eigenvectors are correlated with ruralness? As ruralness became insignificant after incorporation of the spatial eigenvectors we might suspect that it is correlated with a subset of the eigenvectors. Let’s find out which eigenvectors are involved.

kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem %>%
    st_drop_geometry() %>%
    dplyr::select(Laendlichkeit, colnames(sevQp)) %>%
    ggpairs()

Eigenvectors 2 and 1 are negatively correlated with ruralness (“Laendlichkeit”).

It might be that the spatial eigenvectors could be used to moderate the effect of ruralness - i.e. that the regression coefficient varies in space. Likewise we could investigate moderating relationships between other predictors and spatial eigenvectors.

colnames(sevQp)
 [1] "vec2"  "vec6"  "vec4"  "vec10" "vec1"  "vec18" "vec21" "vec11" "vec19"
[10] "vec5" 
modelSevmQpInt <- glm(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
    (vec1 + vec2) * Laendlichkeit + vec4 + vec6 + vec18 + vec21 + vec5 + vec11 +
    vec19 + vec10, family = quasipoisson, offset = log(EWZ), data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem)
summary(modelSevmQpInt)

Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + 
    Hochqualifizierte + (vec1 + vec2) * Laendlichkeit + vec4 + 
    vec6 + vec18 + vec21 + vec5 + vec11 + vec19 + vec10, family = quasipoisson, 
    data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem, 
    offset = log(EWZ))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-56.457  -10.611   -1.923    7.039   64.584  

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -3.7142500  0.0462340 -80.336  < 2e-16 ***
Stimmenanteile.AfD  0.0420443  0.0019723  21.317  < 2e-16 ***
Auslaenderanteil    0.0292675  0.0022308  13.120  < 2e-16 ***
Hochqualifizierte  -0.0048261  0.0018138  -2.661 0.008125 ** 
vec1               -0.2840583  0.1984612  -1.431 0.153158    
vec2                0.8584225  0.2093412   4.101 5.04e-05 ***
Laendlichkeit      -0.0003932  0.0004360  -0.902 0.367675    
vec4               -1.5274509  0.1610968  -9.482  < 2e-16 ***
vec6               -0.8291938  0.1780536  -4.657 4.43e-06 ***
vec18              -0.6106731  0.1490572  -4.097 5.11e-05 ***
vec21              -1.3629828  0.1711811  -7.962 1.94e-14 ***
vec5                0.5642140  0.1649767   3.420 0.000694 ***
vec11              -0.3582028  0.1479799  -2.421 0.015958 *  
vec19              -0.9059447  0.1708088  -5.304 1.92e-07 ***
vec10               0.0975937  0.1615422   0.604 0.546110    
vec1:Laendlichkeit -0.0167823  0.0076383  -2.197 0.028610 *  
vec2:Laendlichkeit  0.0152535  0.0070868   2.152 0.031992 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 221.7727)

    Null deviance: 344479  on 399  degrees of freedom
Residual deviance:  84440  on 383  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4
drop1(modelSevmQpInt, test = "F")
Single term deletions

Model:
sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte + 
    (vec1 + vec2) * Laendlichkeit + vec4 + vec6 + vec18 + vec21 + 
    vec5 + vec11 + vec19 + vec10
                   Df Deviance  F value    Pr(>F)    
<none>                   84440                       
Stimmenanteile.AfD  1   184596 454.2864 < 2.2e-16 ***
Auslaenderanteil    1   121710 169.0505 < 2.2e-16 ***
Hochqualifizierte   1    86012   7.1308 0.0079000 ** 
vec4                1   104417  90.6136 < 2.2e-16 ***
vec6                1    89229  21.7250 4.351e-06 ***
vec18               1    88181  16.9708 4.654e-05 ***
vec21               1    98821  65.2300 8.757e-15 ***
vec5                1    87016  11.6855 0.0006975 ***
vec11               1    85735   5.8741 0.0158275 *  
vec19               1    90660  28.2149 1.845e-07 ***
vec10               1    84521   0.3672 0.5448907    
vec1:Laendlichkeit  1    85511   4.8602 0.0280777 *  
vec2:Laendlichkeit  1    85465   4.6505 0.0316661 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1 - modelSevmQpInt$deviance/modelSevmQpInt$null.deviance
[1] 0.754877
1 - modelSevmQp$deviance/modelSevmQp$null.deviance
[1] 0.7430323

There is a significant interaction between ruralness and spatial eigenvector 2 and 1 while the main effect of ruralness is clearly not significant. However, the explained deviance does not increase much. How does the interaction look like in space?

mRuralVec2 <- kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem %>%
    mutate(vec2Ruralness = vec2 * Laendlichkeit) %>%
    tm_shape() + tm_polygons(col = c("vec2Ruralness"), legend.hist = TRUE, palette = "-RdBu",
    style = "fisher", n = 6, midpoint = 0, legend.reverse = TRUE, title = "Ruralness moderated by eigenvector 2") +
    tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
        legend.outside.size = 0.5, attr.outside = TRUE, legend.outside.position = "left",
        main.title = "Ruralness by sev 2") + tm_scale_bar()

tmap_arrange(mRuralVec2, mRural)

mRuralVec2 <- kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem %>%
    mutate(vec2Ruralness = vec1 * Laendlichkeit) %>%
    tm_shape() + tm_polygons(col = c("vec2Ruralness"), legend.hist = TRUE, palette = "-RdBu",
    style = "fisher", n = 6, midpoint = 0, legend.reverse = TRUE, title = "Ruralness moderated by eigenvector 1") +
    tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
        legend.outside.size = 0.5, attr.outside = TRUE, legend.outside.position = "left",
        main.title = "Ruralness by sev 1") + tm_scale_bar()

tmap_arrange(mRuralVec2, mRural)

lm.morantest(modelSevmQpInt, listw = unionedListw_d)

    Global Moran I for regression residuals

data:  
model: glm(formula = sumCasesTotal ~ Stimmenanteile.AfD +
Auslaenderanteil + Hochqualifizierte + (vec1 + vec2) * Laendlichkeit +
vec4 + vec6 + vec18 + vec21 + vec5 + vec11 + vec19 + vec10, family =
quasipoisson, data =
kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem, offset =
log(EWZ))
weights: unionedListw_d

Moran I statistic standard deviate = 6.9223, p-value = 2.222e-12
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
    1.602309e-01    -2.989182e-06     5.358083e-04 
kreiseCentroids$residSevmQpInt <- residuals(modelSevmQpInt)
kreiseCentroids$residSevmQpIntAbs <- abs(kreiseCentroids$residSevmQpInt)

tm_shape(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded) + tm_polygons(col = "grey") +
    tm_shape(kreiseCentroids) + tm_bubbles(size = "residSevmQpIntAbs", col = "residSevmQpInt",
    palette = "-RdBu", alpha = 0.9, perceptual = TRUE, scale = 0.8, border.alpha = 0.3,
    title.size = "Abs residual", title.col = "Residuals", n = 5) + tm_layout(main.title = "Pearson residuals, SEVM QB\nruralness spatial varying coefficient",
    bg = "darkgrey", legend.outside = TRUE, attr.outside = TRUE) + tm_scale_bar()

modelSevmQpInt2 <- update(modelSevmQpInt, ~. + Auslaenderanteil:vec2)
summary(modelSevmQpInt2)

Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + 
    Hochqualifizierte + vec1 + vec2 + Laendlichkeit + vec4 + 
    vec6 + vec18 + vec21 + vec5 + vec11 + vec19 + vec10 + vec1:Laendlichkeit + 
    vec2:Laendlichkeit + Auslaenderanteil:vec2, family = quasipoisson, 
    data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem, 
    offset = log(EWZ))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-53.822  -10.278   -1.923    7.356   58.824  

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -3.823e+00  4.928e-02 -77.579  < 2e-16 ***
Stimmenanteile.AfD     4.603e-02  2.045e-03  22.512  < 2e-16 ***
Auslaenderanteil       3.303e-02  2.280e-03  14.490  < 2e-16 ***
Hochqualifizierte     -3.396e-03  1.778e-03  -1.910 0.056931 .  
vec1                  -2.644e-01  1.919e-01  -1.377 0.169174    
vec2                   3.858e+00  5.901e-01   6.538 2.00e-10 ***
Laendlichkeit          4.514e-05  4.304e-04   0.105 0.916537    
vec4                  -1.513e+00  1.557e-01  -9.720  < 2e-16 ***
vec6                  -7.843e-01  1.732e-01  -4.529 7.93e-06 ***
vec18                 -6.628e-01  1.443e-01  -4.593 5.94e-06 ***
vec21                 -1.208e+00  1.674e-01  -7.218 2.88e-12 ***
vec5                   5.502e-01  1.597e-01   3.444 0.000636 ***
vec11                 -5.739e-01  1.497e-01  -3.833 0.000148 ***
vec19                 -8.453e-01  1.659e-01  -5.095 5.49e-07 ***
vec10                  3.387e-01  1.628e-01   2.080 0.038173 *  
vec1:Laendlichkeit    -1.555e-02  7.412e-03  -2.097 0.036607 *  
vec2:Laendlichkeit    -1.486e-02  8.818e-03  -1.685 0.092810 .  
Auslaenderanteil:vec2 -2.104e-01  3.875e-02  -5.428 1.02e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 207.0252)

    Null deviance: 344479  on 399  degrees of freedom
Residual deviance:  78251  on 382  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4
1 - modelSevmQpInt2$deviance/modelSevmQpInt2$null.deviance
[1] 0.7728419

If we include the interaction between share of foreigners and sev2 the interaction between ruralness and sev2 becomes insignificant.

modelSevmQpInt2 <- update(modelSevmQpInt2, ~. - vec2:Laendlichkeit)
summary(modelSevmQpInt2)

Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + 
    Hochqualifizierte + vec1 + vec2 + Laendlichkeit + vec4 + 
    vec6 + vec18 + vec21 + vec5 + vec11 + vec19 + vec10 + vec1:Laendlichkeit + 
    Auslaenderanteil:vec2, family = quasipoisson, data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem, 
    offset = log(EWZ))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-55.114  -10.601   -2.052    7.540   60.654  

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -3.8167201  0.0491325 -77.682  < 2e-16 ***
Stimmenanteile.AfD     0.0461689  0.0020456  22.570  < 2e-16 ***
Auslaenderanteil       0.0323683  0.0022466  14.408  < 2e-16 ***
Hochqualifizierte     -0.0035162  0.0017799  -1.975 0.048930 *  
vec1                  -0.2335194  0.1915538  -1.219 0.223564    
vec2                   3.1887116  0.4350967   7.329 1.39e-12 ***
Laendlichkeit          0.0001466  0.0004251   0.345 0.730472    
vec4                  -1.5321164  0.1556024  -9.846  < 2e-16 ***
vec6                  -0.7651248  0.1733036  -4.415 1.32e-05 ***
vec18                 -0.6187969  0.1422214  -4.351 1.74e-05 ***
vec21                 -1.2094141  0.1679333  -7.202 3.17e-12 ***
vec5                   0.5453670  0.1601063   3.406 0.000728 ***
vec11                 -0.5547785  0.1493788  -3.714 0.000234 ***
vec19                 -0.8707075  0.1655418  -5.260 2.40e-07 ***
vec10                  0.2826184  0.1597694   1.769 0.077704 .  
vec1:Laendlichkeit    -0.0116193  0.0070251  -1.654 0.098950 .  
Auslaenderanteil:vec2 -0.1693540  0.0301058  -5.625 3.58e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 207.8777)

    Null deviance: 344479  on 399  degrees of freedom
Residual deviance:  78840  on 383  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4
1 - modelSevmQpInt2$deviance/modelSevmQpInt2$null.deviance
[1] 0.7711321
modelSevmQpInt2 <- update(modelSevmQpInt2, ~. - vec1:Laendlichkeit)
summary(modelSevmQpInt2)

Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + 
    Hochqualifizierte + vec1 + vec2 + Laendlichkeit + vec4 + 
    vec6 + vec18 + vec21 + vec5 + vec11 + vec19 + vec10 + Auslaenderanteil:vec2, 
    family = quasipoisson, data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem, 
    offset = log(EWZ))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-55.008  -10.009   -1.998    7.461   58.275  

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -3.8276546  0.0487983 -78.438  < 2e-16 ***
Stimmenanteile.AfD     0.0469824  0.0019884  23.629  < 2e-16 ***
Auslaenderanteil       0.0327520  0.0022406  14.617  < 2e-16 ***
Hochqualifizierte     -0.0036675  0.0017805  -2.060 0.040093 *  
vec1                  -0.4029965  0.1614259  -2.496 0.012962 *  
vec2                   3.4764206  0.3998051   8.695  < 2e-16 ***
Laendlichkeit          0.0003002  0.0004144   0.724 0.469226    
vec4                  -1.6344266  0.1428255 -11.444  < 2e-16 ***
vec6                  -0.7134649  0.1713180  -4.165 3.86e-05 ***
vec18                 -0.6291195  0.1424453  -4.417 1.31e-05 ***
vec21                 -1.2150457  0.1681284  -7.227 2.69e-12 ***
vec5                   0.5278639  0.1600945   3.297 0.001068 ** 
vec11                 -0.5757764  0.1492522  -3.858 0.000134 ***
vec19                 -0.9210587  0.1630916  -5.647 3.17e-08 ***
vec10                  0.3039040  0.1595634   1.905 0.057579 .  
Auslaenderanteil:vec2 -0.1847810  0.0287029  -6.438 3.62e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 208.6265)

    Null deviance: 344479  on 399  degrees of freedom
Residual deviance:  79408  on 384  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4
1 - modelSevmQpInt2$deviance/modelSevmQpInt2$null.deviance
[1] 0.7694841
modelSevmQpInt2 <- update(modelSevmQpInt2, ~. - Laendlichkeit)
summary(modelSevmQpInt2)

Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + 
    Hochqualifizierte + vec1 + vec2 + vec4 + vec6 + vec18 + vec21 + 
    vec5 + vec11 + vec19 + vec10 + Auslaenderanteil:vec2, family = quasipoisson, 
    data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem, 
    offset = log(EWZ))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-55.793   -9.913   -2.073    7.366   59.329  

Coefficients:
                       Estimate Std. Error  t value Pr(>|t|)    
(Intercept)           -3.805067   0.037489 -101.498  < 2e-16 ***
Stimmenanteile.AfD     0.046717   0.001954   23.908  < 2e-16 ***
Auslaenderanteil       0.032204   0.002108   15.277  < 2e-16 ***
Hochqualifizierte     -0.004154   0.001648   -2.521 0.012113 *  
vec1                  -0.456773   0.143218   -3.189 0.001543 ** 
vec2                   3.384027   0.378769    8.934  < 2e-16 ***
vec4                  -1.639446   0.142565  -11.500  < 2e-16 ***
vec6                  -0.712172   0.171235   -4.159 3.94e-05 ***
vec18                 -0.625980   0.142251   -4.401 1.40e-05 ***
vec21                 -1.222954   0.167732   -7.291 1.76e-12 ***
vec5                   0.525974   0.159953    3.288 0.001101 ** 
vec11                 -0.559127   0.147342   -3.795 0.000172 ***
vec19                 -0.921057   0.163058   -5.649 3.15e-08 ***
vec10                  0.294258   0.158882    1.852 0.064784 .  
Auslaenderanteil:vec2 -0.178247   0.027241   -6.543 1.92e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 208.4964)

    Null deviance: 344479  on 399  degrees of freedom
Residual deviance:  79517  on 385  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4
1 - modelSevmQpInt2$deviance/modelSevmQpInt2$null.deviance
[1] 0.7691665
m1 <- kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem %>%
    mutate(vec2foreigners = vec2 * Auslaenderanteil) %>%
    tm_shape() + tm_polygons(col = c("vec2foreigners"), legend.hist = TRUE, palette = "-RdBu",
    style = "fisher", n = 6, midpoint = 0, legend.reverse = TRUE, title = "Foreigners moderated by eigenvector 2") +
    tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
        legend.outside.size = 0.5, attr.outside = TRUE, legend.outside.position = "left",
        main.title = "Foreigners by sev 2") + tm_scale_bar()

tmap_arrange(m1, mForeigner)

The effect of the share of foreigners differs in space. A potential explanation would be that the effect of the foreign population is moderated by education (at the individual level) as well as by language skills and economic status (of course all factors are interlinked). Also the measures to reach non german speaking inhabitants differs presumably between federal states - leading to potential differences in how those groups obey to rules such as mask wearing, social distancing or with respect to vaccination rates.

Of course there is much more potential to explore relationships between predictors and spatial eigenvectors and we have not even touched on interactions between “normal” predictors or higher order effects. We could (and should) try other predictor eigenvector combinations as well as interactions between the predictors.